.. currentmodule:: brian

.. _example-frompapers-computing with neural synchrony-coincidence detection and synchrony_Fig5E_precision_reliability:

Example: Fig5E_precision_reliability (frompapers/computing with neural synchrony/coincidence detection and synchrony)
=====================================================================================================================

Brette R (2012). Computing with neural synchrony. PLoS Comp Biol. 8(6): e1002561. doi:10.1371/journal.pcbi.1002561
------------------------------------------------------------------------------------------------------------------
Figure 5E. (very long simulation)

Caption (Fig 5E). Precision and reliability of spike timing as a function of SNR.

Simulations are run in parallel on all cores but one.

::

    from brian import *
    import multiprocessing
    
    def autocor(PSTH,N=None,T=20*ms,bin=None):
        '''
        Autocorrelogram of PSTH, to calculate a shuffled autocorrelogram
        
        N = number of spike trains
        T = temporal window
        bin = PSTH bin
        The baseline is not subtracted.
        
        Returns times,SAC
        '''
        if bin is None:
            bin=defaultclock.dt
        n=len(PSTH)
        p=int(T/bin)
        SAC=zeros(p)
        if N is None:
            SAC[0]=mean(PSTH*PSTH)
        else: # correction to exclude self-coincidences
            PSTHnoself=clip(PSTH-1./(bin*N),0,Inf)
            SAC[0]=mean(PSTH*PSTHnoself)*N/(N-1.)
        SAC[1:]=[mean(PSTH[:-i]*PSTH[i:]) for i in range(1,p)]
        SAC=hstack((SAC[::-1],SAC[1:]))
        return (arange(len(SAC))-len(SAC)/2)*bin,SAC
    
    def halfwidth(x):
        '''
        Returns half-width of function given by x in bin numbers.
        This is used to calculate the precision (left panel).
        '''
        M,n=max(x),argmax(x)
        return find(x[n:]<M/2)[0]+n-find(x[:n]<M/2)[-1]
    
    def reproducibility(SNR):
        '''
        Calculates the precision (timescale) and reliability (strength) for a given
        signal-to-noise ratio.
        '''
        sys.stdout.write("SNR:"+str(SNR)+'\n')
        sys.stdout.flush() # we use this instead of print because of multiprocessing
        reinit_default_clock() # important because we do multiple simulations
        # The common noisy input
        N=5000                 # number of neurons simultaneously simulated
        duration=30*second     # duration of one simulation, 200 seconds in the paper
        tau_noise=5*ms
        input=NeuronGroup(1,model='dx/dt=-x/tau_noise+(2./tau_noise)**.5*xi:1')
        
        # The noisy neurons receiving the same input
        tau=10*ms
        sigma=.5 # input amplitude
        Z=sigma*sqrt((tau_noise+tau)/(tau_noise*(SNR**2+1))) # normalizing factor
        eqs_neurons='''
        dx/dt=(Z*(SNR*I+u)-x)/tau:1
        du/dt=-u/tau_noise+(2./tau_noise)**.5*xi:1
        I : 1
        '''
        neurons=NeuronGroup(N,model=eqs_neurons,threshold=1,reset=0,refractory=5*ms)
        neurons.x=rand(N) # random initial conditions
        neurons.I=linked_var(input,'x')
        rate=PopulationRateMonitor(neurons) # PSTH
        
        run(duration)
        
        t,SAC=autocor(rate.rate,N,T=30*ms)
        timescale=float(halfwidth(SAC-mean(rate.rate)**2))*defaultclock.dt # precision
        strength=sum(SAC-mean(rate.rate)**2)*float(defaultclock.dt)/mean(rate.rate) # reliability
        
        return timescale,strength
    
    if __name__=='__main__':
        pool = multiprocessing.Pool(multiprocessing.cpu_count()-1) # all cores but one
        SNRdB= linspace(-10,15,20) # 100 points in the paper
        SNR = 10.**(SNRdB/10.)
        results = pool.map(reproducibility, SNR) # launches multiple processes
        timescale,strength=zip(*results)
    
        # Figure
        subplot(211)
        plot(SNRdB,timescale*1000)
        xlabel('SNR (dB)')
        ylabel('Precision (ms)')
        subplot(212)
        plot(SNRdB,strength*100)
        xlabel('SNR (dB)')
        ylabel('Reliability (%)')
        show()